function [mm_D1] = DO_Appendix_D1(param)
%% ==================== Set Parameters and Grids ========================%

    %----- Specify Parameters -----%
    nu = param(1);
    beta = param(2);            % Discount factor (factors in exit probability)

    alpha = param(3);           % Production function parameter
    d = param(4);               % Human capital depreciation rate
    ell = param(5);             % Matching function parameter
    lambda = param(6);          % Probability of search OTJ
    chi = param(7);             % Investment cost parameter
    b=param(8);                 % Unemployment benefit

    zetaF =param(9);              % --- Equivalent to zeta^F in Appendix D
    c =param(10);               % Cost of posting a vacancy
    delta =param(11);           % Separation rate
    h_sigma =param(12);         % Standard deviation of initial human capital distribution

    pr_phi(1) = param(13);      % Probability of entering the model with phi=1 (low fixed productivity type)
    pr_phi(2) = param(14);      % Probability of entering the model with phi=2 (high fixed productivity type)

    z(1) = param(15);           % Productivity parameter f(phi,h) = z(phi)h^(alpha)
    z(2) = param(16);
    gamma=2;
    vfi_toler = 0.0001;
    dist_toler = 0.0000001;
    
    pr_phi(1) = 1.0 - 0.3084251;  % Probability of entering the model with phi=1 (low fixed productivity type)
    pr_phi(2) = 0.3084251;      % Probability of entering the model with phi=2 (high fixed productivity type)
    
    z(1) = 1;                 % Productivity parameter f(phi,h) = z(phi)h^(alpha)
    z(2) = 1.559533;
    
    h_mu = 1;                 % Mean of initial human capital distribution
    nphi = 2;                 % Number of possible phi values (fixed productivity types)
    
    
    hlb = 0.0;               % Human capital grid lower bound 
    hub = 7.0;               % Human capital grid upper bound 
    nh = 3000;                % Number of possible human capital values
    h_grid = linspace(hlb,hub,nh); % Create human capital grid
    
    % ======================= Create discretized h distribution ==================== %
    cdf_h = zeros(1,nh);      % CDF of initial h distribution
    pdf_h = zeros(1,nh);      % PDF of initial h distribution
    
    h_bin_size = (hub-hlb)/(nh-1);
    for ind_h = 1:nh
        if (ind_h ==1)
            bin_ub_now = h_grid(ind_h) + (h_bin_size)/2;
            cdf_h(ind_h) = normcdf(bin_ub_now,h_mu,h_sigma);
            pdf_h(ind_h) = normcdf(bin_ub_now,h_mu,h_sigma);
        else
            bin_ub_now = h_grid(ind_h) + (h_bin_size)/2;
            bin_ub_last = h_grid(ind_h-1) + (h_bin_size)/2;
            cdf_h(ind_h) = normcdf(bin_ub_now,h_mu,h_sigma);
            pdf_h(ind_h) = normcdf(bin_ub_now,h_mu,h_sigma) - normcdf(bin_ub_last,h_mu,h_sigma);
        end
    end
    %figure
    %plot(h_grid,pdf_h)
    
    
    %% ============== Backward Loop over Ages to Get Value and Policy Functions ================ %%
    %---- Prespecify Value Functions
    
    U = zeros(nphi,nh);                % Value of unemployment in production phase                     
    U_hat  = zeros(nphi,nh);           % Value of unemployment in search phase - when searching for optimal job type
    U_hat_update  = zeros(nphi,nh);
    
    VN = zeros(nphi,nh);               % Value of employment in non-outsourced (N) job in production phase
    VO = zeros(nphi,nh);               % Value of employment in outsourced (O) job in production phase
    VN_hat = zeros(nphi,nh);           % Value of employment in N job in search phase
    VO_hat = zeros(nphi,nh);           % Value of employment in O job in search phase
    VN_hat_update = zeros(nphi,nh);
    VO_hat_update = zeros(nphi,nh);
    
    %---- Prespecify Value Functions
    opt_ind_hp_N = ones(nphi,nh);     % Optimal index for next period h value (h') if in non-outsourced job
    opt_ind_hp_O = ones(nphi,nh);     % Opitmal index for next period h value (h') if in outsourced job
    
    GU = zeros(nphi,nh);              % Indicates optimal job type to search for if unemployed - 1=non-outsourced (N) and 2=outsourced (O)
    GN = zeros(nphi,nh);              % Indicates optimal job type to search for if employed in N job - 1=non-outsourced and 2=outsourced
    GO = zeros(nphi,nh);              % Indicates optimal job type to search for if employed in O job - 1=non-outsourced and 2=outsourced
    
    opt_p_U = zeros(nphi,nh);         % Optimal job-finding rate associated with optimal search choice when unemployed
    opt_p_N = zeros(nphi,nh);         % Optimal job-finding rate associated with optimal search choice when employed in N job
    opt_p_O = zeros(nphi,nh);         % Optimal job-finding rate associated with optimal search choice when employed in O job
    
    opt_x_U = zeros(nphi,nh);         % x contract value associated with optimal search choice when unemployed
    opt_x_N = zeros(nphi,nh);         % x contract value associated with optimal search choice when employed in N job
    opt_x_O = zeros(nphi,nh);         % x contract value associated with optimal search choice when employed in O job
    
    opt_eta_U = zeros(nphi,nh);       % eta (surplus rate) associated with optimal search choice when unemployed
    opt_eta_N = zeros(nphi,nh);       % eta (surplus rate) associated with optimal search choice when employed in N job 
    opt_eta_O = zeros(nphi,nh);       % eta (surplus rate) associated with optimal search choice when employed in O job
    
    opt_theta_U = zeros(nphi,nh);     % theta (market tightness) associated with optimal search choice when unemployed
    opt_theta_N = zeros(nphi,nh);     % theta (market tightness) associated with optimal search choice when employed in N job 
    opt_theta_O = zeros(nphi,nh);     % theta (market tightness) associated with optimal search choice when employed in O job
    
    % Find index for next period human capital (h') if the agent only lets human capital depreciate
    ind_hp_min = ones(1,nh);          % Minimum index for h' (index for h' if human capital only depreciates)
    for ind_h = 1:nh
        h = h_grid(ind_h);
        error_min = 100;
        for ind_hp = 1:nh
            hp = h_grid(ind_hp);       % hp is the value of h'
            error = abs(hp - h*(1-d)); % Find index that minimizes (hp-h*(1-d))
            if (error<error_min)
                error_min = error;
                ind_hp_min(ind_h) = ind_hp;
            end
        end
    end
    
    %--------------- Begin Value Function Iteration ----------------%
    vfi_error = 1; % Set vfi_error > vfi_toler
    while (vfi_error > vfi_toler)
    
        for ind_phi = 1:nphi          % Loop over phi types
            for ind_h = 1:nh          % Loop over current h values
                h = h_grid(ind_h);    % h = current human capital (h) value 
                
                b_effect = b*z(ind_phi)*(h^alpha);  % Unemployment value
                c_effect = c;                                % Search cost
    
                % ---------- Solve hc Investment Decision ----------%
                max_VN = -9999;
                max_VO = -9999;
    
                for ind_hp = ind_hp_min(ind_h):nh      % Loop over possible next period hp (h') values
                    hp = h_grid(ind_hp);
    
                    temp_VN = z(ind_phi)*(h^alpha) - ((hp - h*(1-d))*chi)^gamma ...
                        + beta*(1-delta)*VN_hat(ind_phi,ind_hp) + beta*delta*U_hat(ind_phi,ind_hp);
    
                    temp_VO = (1-zetaF)*(z(ind_phi)*(h^alpha)) - ((hp - h*(1-d))*chi)^gamma ...
                        + beta*(1-delta)*VO_hat(ind_phi,ind_hp) + beta*delta*U_hat(ind_phi,ind_hp);
    
                    if (temp_VN>max_VN)                % Save ind_hp where N job employment value is highest
                        max_VN = temp_VN;
                        VN(ind_phi,ind_h) = temp_VN;
                        opt_ind_hp_N(ind_phi,ind_h) = ind_hp;
                    end
                    if (temp_VO>max_VO)                % Save ind_hp where O job employment value is highest
                        max_VO = temp_VO;
                        VO(ind_phi,ind_h) = temp_VO;
                        opt_ind_hp_O(ind_phi,ind_h) = ind_hp;
                    end
    
                end % end loop over hp options
                U(ind_phi,ind_h) = b_effect + beta*U_hat(ind_phi,ind_hp_min(ind_h));
    
                %---------------------------------------------------%
                % ------- Solve Search Decision of Unemployed ----- %
                xO = VO(ind_phi,ind_h) - c_effect;
                etaO = (xO - U(ind_phi,ind_h))/(VO(ind_phi,ind_h) - U(ind_phi,ind_h));
                t_Uhat_O = xO;
    
                thetaN = (((VN(ind_phi,ind_h) - U(ind_phi,ind_h))/c_effect)^(ell/(1+ell)) - 1)^(1/ell);
                if (~isreal(thetaN))
                    pN = 0;
                    xN = 0;
                    etaN = 0;
                else
                    pN = thetaN/((1+thetaN^ell)^(1/ell));
                    qN = pN/thetaN;
                    xN = VN(ind_phi,ind_h) - (c_effect/qN);
                    etaN = (xN - U(ind_phi,ind_h))/(VN(ind_phi,ind_h) - U(ind_phi,ind_h));
                end
                t_Uhat_N = pN*xN + (1-pN)*U(ind_phi,ind_h);
    
                % Now Consider which job type's value is highest
                if (t_Uhat_N>t_Uhat_O)
                    GU(ind_phi,ind_h) = 1;
                    opt_p_U(ind_phi,ind_h) = pN;
                    U_hat_update(ind_phi,ind_h) = t_Uhat_N;
                    opt_x_U(ind_phi,ind_h) = xN;
                    opt_eta_U(ind_phi,ind_h) = etaN;
                    opt_theta_U(ind_phi,ind_h) = thetaN;
                elseif (t_Uhat_O>t_Uhat_N)
                    GU(ind_phi,ind_h) = 2;
                    opt_p_U(ind_phi,ind_h) = 1;
                    U_hat_update(ind_phi,ind_h) = t_Uhat_O;
                    opt_x_U(ind_phi,ind_h) = xO;
                    opt_eta_U(ind_phi,ind_h) = etaO;
                    opt_theta_U(ind_phi,ind_h) = 1.0;
                else
                    GU(ind_phi,ind_h) = 3;
                    opt_p_U(ind_phi,ind_h) = pN;
                    U_hat_update(ind_phi,ind_h) = t_Uhat_N;
                    opt_x_U(ind_phi,ind_h) = xN;
                    opt_eta_U(ind_phi,ind_h) = etaN;
                    opt_theta_U(ind_phi,ind_h) = thetaN;
                end
    
                %---------------------------------------------------%
                % ------- Solve Search Decision of Employed --------%
                xO = VO(ind_phi,ind_h) - c_effect;
                etaO = (xO - U(ind_phi,ind_h))/(VO(ind_phi,ind_h) - U(ind_phi,ind_h));
    
                t_VNhat_O = lambda*xO + (1-lambda)*VN(ind_phi,ind_h);
                t_VOhat_O = lambda*xO + (1-lambda)*VO(ind_phi,ind_h);
    
                thetaN_N = (((VN(ind_phi,ind_h) - VN(ind_phi,ind_h))/c_effect)^(ell/(1+ell)) - 1)^(1/ell);
                if (~isreal(thetaN_N))
                    pN_N = 0;
                    xN_N = 0;
                    etaN_N = 0;
                else
                    pN_N = thetaN_N/((1+thetaN_N^ell)^(1/ell));
                    qN_N = pN_N/thetaN_N;
                    xN_N = VN(ind_phi,ind_h) - (c_effect/qN_N);
                    etaN_N = (xN_N - U(ind_phi,ind_h))/(VN(ind_phi,ind_h) - U(ind_phi,ind_h));
                end
                t_VNhat_N = lambda*pN_N*xN_N + (1-lambda*pN_N)*VN(ind_phi,ind_h);
    
                thetaO_N = (((VN(ind_phi,ind_h) - VO(ind_phi,ind_h))/c_effect)^(ell/(1+ell)) - 1)^(1/ell);
                if (~isreal(thetaO_N))
                    pO_N = 0;
                    xO_N = 0;
                    etaO_N = 0;
                else
                    pO_N = thetaO_N/((1+thetaO_N^ell)^(1/ell));
                    qO_N = pO_N/thetaO_N;
                    xO_N = VN(ind_phi,ind_h) - (c_effect/qO_N);
                    etaO_N = (xO_N - U(ind_phi,ind_h))/(VN(ind_phi,ind_h) - U(ind_phi,ind_h));
                end  
                t_VOhat_N = lambda*pO_N*xO_N + (1-lambda*pO_N)*VO(ind_phi,ind_h);
    
                % For those Currently in N - Now Consider which job type is highest
                if (t_VNhat_N>t_VNhat_O)
                    GN(ind_phi,ind_h) = 1;
                    opt_p_N(ind_phi,ind_h) = pN_N;
                    VN_hat_update(ind_phi,ind_h) = t_VNhat_N;
                    opt_x_N(ind_phi,ind_h) = xN_N;
                    opt_eta_N(ind_phi,ind_h) = etaN_N;
                    opt_theta_N(ind_phi,ind_h) = thetaN_N;
                elseif (t_VNhat_O>t_VNhat_N)
                    GN(ind_phi,ind_h) = 2;
                    opt_p_N(ind_phi,ind_h) = 1;
                    VN_hat_update(ind_phi,ind_h) = t_VNhat_O;
                    opt_x_N(ind_phi,ind_h) = xO;
                    opt_eta_N(ind_phi,ind_h) = etaO;
                    opt_theta_N(ind_phi,ind_h) = 1.0;
                else
                    GN(ind_phi,ind_h) = 3;
                    opt_p_N(ind_phi,ind_h) = pN_N;
                    VN_hat_update(ind_phi,ind_h) = t_VNhat_N;
                    opt_x_N(ind_phi,ind_h) = xN_N;
                    opt_eta_N(ind_phi,ind_h) = etaN_N;
                    opt_theta_N(ind_phi,ind_h) = thetaN_N;
                end
    
                % For those Currently in O - Now Consider which job type is highest
                if (t_VOhat_N>t_VOhat_O)
                    GO(ind_phi,ind_h) = 1;
                    opt_p_O(ind_phi,ind_h) = pO_N;
                    VO_hat_update(ind_phi,ind_h) = t_VOhat_N;
                    opt_x_O(ind_phi,ind_h) = xO_N;
                    opt_eta_O(ind_phi,ind_h) = etaO_N;
                    opt_theta_O(ind_phi,ind_h) = thetaO_N;
                elseif (t_VOhat_O>t_VOhat_N)
                    GO(ind_phi,ind_h) = 2;
                    opt_p_O(ind_phi,ind_h) = 1;
                    VO_hat_update(ind_phi,ind_h) = t_VOhat_O;
                    opt_x_O(ind_phi,ind_h) = xO;
                    opt_eta_O(ind_phi,ind_h) = etaO;
                    opt_theta_O(ind_phi,ind_h) = 1.0;
                else 
                    GO(ind_phi,ind_h) = 3;
                    opt_p_O(ind_phi,ind_h) = pO_N;
                    VO_hat_update(ind_phi,ind_h) = t_VOhat_N;
                    opt_x_O(ind_phi,ind_h) = xO_N;
                    opt_eta_O(ind_phi,ind_h) = etaO_N;
                    opt_theta_O(ind_phi,ind_h) = thetaO_N;
                end
    
            end % End loop over h values
        end % End loop over phi values
    
        error_U_hat = max(max(abs(U_hat_update(:,:) - U_hat(:,:))))/max(max(abs(U_hat_update(:,:))));
        error_VN_hat = max(max(abs(VN_hat_update(:,:) - VN_hat(:,:))))/max(max(abs(VN_hat_update(:,:))));
        error_VO_hat = max(max(abs(VO_hat_update(:,:) - VO_hat(:,:))))/max(max(abs(VO_hat_update(:,:))));
    
        vfi_error = max([error_U_hat,error_VN_hat,error_VO_hat])
    
        U_hat(:,:) = U_hat_update(:,:);
        VN_hat(:,:) = VN_hat_update(:,:);
        VO_hat(:,:) = VO_hat_update(:,:);
        
    end % End while (vfi_error > vfi_toler)  
    
    %==================================================================================================================%
    %% ================================ Solve for Steady State Distribution ========================================= %%
    %==================================================================================================================%
    
    dist_U = zeros(nphi,nh);           % Steady state mass of unemployed agents
    dist_N = zeros(nphi,nh);           % Steady state mass of agents employed in non-outsourced job
    dist_O = zeros(nphi,nh);           % Steady state mass of agents employed in outsourced job
    
    dist_U_update = zeros(nphi,nh);
    dist_N_update = zeros(nphi,nh);
    dist_O_update = zeros(nphi,nh);
    
    dist_U(1,1) = pr_phi(1);           % Start with a guess regarding how the mass of agents is distributed
    dist_U(2,1) = pr_phi(2);
    
    dist_error = 1;
    while (dist_error>dist_toler)      % Iterate until distribution guess equals resulting distribution
    
        for ind_phi = 1:nphi           % Loop over phi types
            for ind_h = 1:nh           % Loop over current h values 
                %=================================%
                %========= New Entrants ==========%
                %=================================%
                % New Entrant - Finds Job
                if (GU(ind_phi,ind_h)==1 || GU(ind_phi,ind_h)==3)     % Searches for N
                    dist_N_update(ind_phi,ind_h) = dist_N_update(ind_phi,ind_h)...
                        + opt_p_U(ind_phi,ind_h)*pr_phi(ind_phi)*nu*pdf_h(ind_h);
                else                                                   % Searches for O
                    dist_O_update(ind_phi,ind_h) = dist_O_update(ind_phi,ind_h)...
                        + opt_p_U(ind_phi,ind_h)*pr_phi(ind_phi)*nu*pdf_h(ind_h);
                end
                % New Entrant - Does Not Find Job
                dist_U_update(ind_phi,ind_h) = dist_U_update(ind_phi,ind_h) + (1-opt_p_U(ind_phi,ind_h))*pr_phi(ind_phi)*nu*pdf_h(ind_h);
    
                %=================================%
                %===== Previously Unemployed =====%
                %=================================%
                ind_hp = ind_hp_min(ind_h);             % New h value for the unemployed is what they get from allowing human capital to depreciate
    
                % Unemployed - Finds a job
                if (GU(ind_phi,ind_hp) == 1 || GU(ind_phi,ind_hp) == 3)     % Searches for N
                    dist_N_update(ind_phi,ind_hp) = dist_N_update(ind_phi,ind_hp)...
                        + opt_p_U(ind_phi,ind_hp)*dist_U(ind_phi,ind_h)*(1-nu);
                else                                                         % Searches for O
                    dist_O_update(ind_phi,ind_hp) = dist_O_update(ind_phi,ind_hp)...
                        + opt_p_U(ind_phi,ind_hp)*dist_U(ind_phi,ind_h)*(1-nu);
                end 
                % Unemployed - Does not find a job
                dist_U_update(ind_phi,ind_hp) = dist_U_update(ind_phi,ind_hp) + (1-opt_p_U(ind_phi,ind_hp))*dist_U(ind_phi,ind_h)*(1-nu);
    
                %=================================%
                %======= Previously in N =========%
                %=================================%
                ind_hp = opt_ind_hp_N(ind_phi,ind_h);   % New h value for employed in N job is the result of their optimal investment choice
    
                % ----- N - No Separation Shock - Finds New Job
                if (GN(ind_phi,ind_hp)==1 || GN(ind_phi,ind_hp)==3)      % N - Search OTJ for N
                    dist_N_update(ind_phi,ind_hp) = dist_N_update(ind_phi,ind_hp)...
                        + lambda*opt_p_N(ind_phi,ind_hp)*(1-delta)*dist_N(ind_phi,ind_h)*(1-nu);
                else                                                      % N - Search OTJ for O
                    dist_O_update(ind_phi,ind_hp) = dist_O_update(ind_phi,ind_hp)...
                        + lambda*opt_p_N(ind_phi,ind_hp)*(1-delta)*dist_N(ind_phi,ind_h)*(1-nu);
                end
    
                % ----- N - No Separation Shock - Does Not Find New Job
                dist_N_update(ind_phi,ind_hp) = dist_N_update(ind_phi,ind_hp)...
                    + (1-lambda*opt_p_N(ind_phi,ind_hp))*(1-delta)*dist_N(ind_phi,ind_h)*(1-nu);
    
                % ----- N - Separation Shock - Finds Job
                if (GU(ind_phi,ind_hp)==1 || GU(ind_phi,ind_hp)==3)       % N - Separation - Search for N
                    dist_N_update(ind_phi,ind_hp) = dist_N_update(ind_phi,ind_hp)...
                        + opt_p_U(ind_phi,ind_hp)*delta*dist_N(ind_phi,ind_h)*(1-nu);
                else                                                       % N - Separation - Search for O
                    dist_O_update(ind_phi,ind_hp) = dist_O_update(ind_phi,ind_hp)...
                        + opt_p_U(ind_phi,ind_hp)*delta*dist_N(ind_phi,ind_h)*(1-nu);
                end
                
                % ----- N - Separation Shock - Does Not Find New Job
                dist_U_update(ind_phi,ind_hp) = dist_U_update(ind_phi,ind_hp)...
                    + (1-opt_p_U(ind_phi,ind_hp))*delta*dist_N(ind_phi,ind_h)*(1-nu);
    
                %=================================%
                %======= Previously in O =========%
                %=================================%
                ind_hp = opt_ind_hp_O(ind_phi,ind_h);   % New h value for employed in O job is the result of their optimal investment choice
    
                % ----- O - No Separation Shock - Finds New Job
                if (GO(ind_phi,ind_hp)==1 || GO(ind_phi,ind_hp)==3)      % O - Search OTJ for N
                    dist_N_update(ind_phi,ind_hp) = dist_N_update(ind_phi,ind_hp)...
                        + lambda*opt_p_O(ind_phi,ind_hp)*(1-delta)*dist_O(ind_phi,ind_h)*(1-nu);
                else                                                      % O - Search OTJ for O
                    dist_O_update(ind_phi,ind_hp) = dist_O_update(ind_phi,ind_hp)...
                        + lambda*opt_p_O(ind_phi,ind_hp)*(1-delta)*dist_O(ind_phi,ind_h)*(1-nu);
                end
                % ----- O - No Separation Shock - Does Not Find New Job
                dist_O_update(ind_phi,ind_hp) = dist_O_update(ind_phi,ind_hp)...
                    + (1-lambda*opt_p_O(ind_phi,ind_hp))*(1-delta)*dist_O(ind_phi,ind_h)*(1-nu);
    
                % ----- O - Separation Shock - Finds Job
                if (GU(ind_phi,ind_hp)==1 || GU(ind_phi,ind_hp)==3)       % O - Separation - Search for N
                    dist_N_update(ind_phi,ind_hp) = dist_N_update(ind_phi,ind_hp)...
                        + opt_p_U(ind_phi,ind_hp)*delta*dist_O(ind_phi,ind_h)*(1-nu);
                else                                                       % O - Separation - Search for O
                    dist_O_update(ind_phi,ind_hp) = dist_O_update(ind_phi,ind_hp)...
                        + opt_p_U(ind_phi,ind_hp)*delta*dist_O(ind_phi,ind_h)*(1-nu);
                end
    
                % ----- O - Separation Shock - Does Not Find New Job
                dist_U_update(ind_phi,ind_hp) = dist_U_update(ind_phi,ind_hp)...
                    + (1-opt_p_U(ind_phi,ind_hp))*delta*dist_O(ind_phi,ind_h)*(1-nu);
    
            end % End loop over last period h value
        end % End loop over phi values
    
        error_dist_U = max(max(abs(dist_U_update(:,:) - dist_U(:,:))));
        error_dist_N = max(max(abs(dist_N_update(:,:) - dist_N(:,:))));
        error_dist_O = max(max(abs(dist_O_update(:,:) - dist_O(:,:))));
    
        dist_error = max([error_dist_U,error_dist_N,error_dist_O]);
    
        dist_U(:,:) = dist_U_update(:,:);
        dist_N(:,:) = dist_N_update(:,:);
        dist_O(:,:) = dist_O_update(:,:);
    
        dist_U_update(:,:) = 0;
        dist_N_update(:,:) = 0;
        dist_O_update(:,:) = 0;
    
    
    end % end while (dist_error>dist_toler)
    

    %% =============================================================================================================== %%
    %======================================== Simulation for Wage Moments and Results ==================================%
    %================================================================================================================= %%
    
    % Need to follow "respondents" from entrance because SS distribution does not save the eta values %
    simul_size = 20000;     % Number of individuals to follow in simulated sample
    nt = 400;               % Number of quarters to follow each individual - 100 yrs max - most will exit the model before this
    simul(simul_size,nt) = struct;
    rng(123456789);         % Set random number seed so results are exactly the same each time the code is ran   
    
    for id = 1:simul_size   % Loop over each individual in the sample
        id
        for t = 1:nt        % Loop over each time period the individual is observed
            
            if (t==1)
                % For new entrants, assign phi and h, then see if they are successful in their job search from unemployment %
                % Assign starting h value
                r1 = rand;
                simul(id,t).h_ind = 1;  % If r1 is not bigger than any CDF values, assign to the lowest slot
                for ind_h = 2:nh
                    if (r1>cdf_h(ind_h-1)) 
                        simul(id,t).h_ind = ind_h;
                    end 
                end
                ind_h = simul(id,t).h_ind;
                h = h_grid(ind_h);
                % Now assign phi type
                r1 = rand;
                if r1 < pr_phi(1)
                    simul(id,t).phi_ind = 1;
                else
                    simul(id,t).phi_ind = 2;
                end
                ind_phi = simul(id,t).phi_ind;
                % See if new entrant is successful in job search from unemployment - - If successful, save their wage %
                r1 = rand;
                if (r1 < opt_p_U(ind_phi,ind_h)) % Finds job
                    if (GU(ind_phi,ind_h)==1 || GU(ind_phi,ind_h)==3)        % Search for N
                        simul(id,t).stat = 1;                                % Status "stat" 1 = in non-outsourced job
                        simul(id,t).eta = opt_eta_U(ind_phi,ind_h);
                        eta = simul(id,t).eta;
                        
                        ind_hp = opt_ind_hp_N(ind_phi,ind_h);                % This is the h' value they will choose for next period 
                        hp = h_grid(ind_hp);
                        simul(id,t).wage = eta*(VN(ind_phi,ind_h) - U(ind_phi,ind_h)) + U(ind_phi,ind_h) + (chi*(hp - h*(1-d)))^gamma...
                            - beta*delta*U_hat(ind_phi,ind_hp) - beta*(1-delta)*lambda*opt_p_N(ind_phi,ind_hp)*opt_x_N(ind_phi,ind_hp)...
                            - beta*(1-delta)*(1-lambda*opt_p_N(ind_phi,ind_hp))*(eta*VN(ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp));
                            
                    else                                                      % Search for O
                        simul(id,t).stat = 2;                                 % Status "stat" 2 = in outsourced job 
                        simul(id,t).eta = opt_eta_U(ind_phi,ind_h);
                        eta = simul(id,t).eta;
                        
                        ind_hp = opt_ind_hp_O(ind_phi,ind_h);                % This is the h' value they will choose for next period 
                        hp = h_grid(ind_hp);
     
                        simul(id,t).wage = eta*(VO(ind_phi,ind_h) - U(ind_phi,ind_h)) + U(ind_phi,ind_h) + (chi*(hp - h*(1-d)))^gamma...
                            - beta*delta*U_hat(ind_phi,ind_hp) - beta*(1-delta)*lambda*opt_p_O(ind_phi,ind_hp)*opt_x_O(ind_phi,ind_hp)...
                            - beta*(1-delta)*(1-lambda*opt_p_O(ind_phi,ind_hp))*(eta*VO(ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp));   
                    end   
                else  % Does not find job
                    simul(id,t).stat = 0;                                    % Status "stat" 0 = unemployed  
                    simul(id,t).wage = 0.0;
                    simul(id,t).eta = 0.0;
                end   % End if t=1 individual finds job or not
                
            else % so t>1
                simul(id,t).phi_ind = simul(id,t-1).phi_ind;                 % Individual's phi value is permanent
                ind_phi = simul(id,t).phi_ind;
                
                if (simul(id,t-1).stat==999) % if dead                       % If the agent already exited the model they remain dead             
                    simul(id,t).stat=999;                                    % Status "stat" 999 = dead
                else   % was not dead in last period
                    rnu = rand;
                    if (rnu<nu) % Dies at end of last period
                        simul(id,t).stat=999;
                    else        % Does NOT die at end of last period
                        
                        if ((simul(id,t-1).stat)==0)        % If unemployed last period 
                            % Recall h from last period and apply any depreciation
                            ind_h = simul(id,t-1).h_ind;
                            simul(id,t).h_ind =  ind_hp_min(ind_h);
                            ind_h = simul(id,t).h_ind;
                            h = h_grid(ind_h);
                            
                            % If unemployed last period - see if matches from unemployment !
                            r1 = rand;
                            if (r1 < opt_p_U(ind_phi,ind_h))    % Finds job
                                if (GU(ind_phi,ind_h)==1 || GU(ind_phi,ind_h)==3)      % Search for N job
                                    simul(id,t).stat = 1; 
                                    simul(id,t).eta = opt_eta_U(ind_phi,ind_h);
                                    eta = simul(id,t).eta;
                                    ind_hp = opt_ind_hp_N(ind_phi,ind_h);
                                    hp = h_grid(ind_hp);
                                    simul(id,t).wage = eta*(VN(ind_phi,ind_h) - U(ind_phi,ind_h)) + U(ind_phi,ind_h) + (chi*(hp - h*(1-d)))^gamma...
                                        - beta*delta*U_hat(ind_phi,ind_hp) - beta*(1-delta)*lambda*opt_p_N(ind_phi,ind_hp)*opt_x_N(ind_phi,ind_hp)...
                                        - beta*(1-delta)*(1-lambda*opt_p_N(ind_phi,ind_hp))*(eta*VN(ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp));  
                                else                                                    % Search for O job
                                    simul(id,t).stat = 2; 
                                    simul(id,t).eta = opt_eta_U(ind_phi,ind_h);
                                    eta = simul(id,t).eta;
                                    ind_hp = opt_ind_hp_O(ind_phi,ind_h);
                                    hp = h_grid(ind_hp);
                                    simul(id,t).wage = eta*(VO(ind_phi,ind_h) - U(ind_phi,ind_h)) + U(ind_phi,ind_h) + (chi*(hp - h*(1-d)))^gamma...
                                        - beta*delta*U_hat(ind_phi,ind_hp) - beta*(1-delta)*lambda*opt_p_O(ind_phi,ind_hp)*opt_x_O(ind_phi,ind_hp)...
                                        - beta*(1-delta)*(1-lambda*opt_p_O(ind_phi,ind_hp))*(eta*VO(ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp));  
                                end
                            else                                 % Unemployed last period and does not find job%
                                simul(id,t).stat = 0;
                                simul(id,t).wage = 0.0;
                                simul(id,t).eta = 0.0;
                            end
     
                        elseif ((simul(id,t-1).stat)==1)    % If in N last period  
                            % Recall h from last period and apply any investment
                            ind_h = simul(id,t-1).h_ind;
                            simul(id,t).h_ind = opt_ind_hp_N(ind_phi,ind_h);
                            ind_h = simul(id,t).h_ind;
                            h = h_grid(ind_h);
                            
                            % If N last period - see if separation shock %
                            rdelta = rand;
                            if (rdelta < delta)  % N last period - separation shock    %
                                % See if finds job from unemployment
                                r1 = rand;
                                if (r1 < opt_p_U(ind_phi,ind_h)) % N last period - separation shock - finds job from unemployment
                                    if (GU(ind_phi,ind_h)==1 || GU(ind_phi,ind_h)==3)      % Search for N
                                        simul(id,t).stat = 1; 
                                        simul(id,t).eta = opt_eta_U(ind_phi,ind_h);
                                        eta = simul(id,t).eta;
                                        ind_hp = opt_ind_hp_N(ind_phi,ind_h);
                                        hp = h_grid(ind_hp);
                                        simul(id,t).wage = eta*(VN(ind_phi,ind_h) - U(ind_phi,ind_h)) + U(ind_phi,ind_h) + (chi*(hp - h*(1-d)))^gamma...
                                            - beta*delta*U_hat(ind_phi,ind_hp) - beta*(1-delta)*lambda*opt_p_N(ind_phi,ind_hp)*opt_x_N(ind_phi,ind_hp)...
                                            - beta*(1-delta)*(1-lambda*opt_p_N(ind_phi,ind_hp))*(eta*VN(ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp));  
                                    else                                                    % Search for O
                                        simul(id,t).stat = 2; 
                                        simul(id,t).eta = opt_eta_U(ind_phi,ind_h);
                                        eta = simul(id,t).eta;
                                        ind_hp = opt_ind_hp_O(ind_phi,ind_h);
                                        hp = h_grid(ind_hp);
                                        simul(id,t).wage = eta*(VO(ind_phi,ind_h) - U(ind_phi,ind_h)) + U(ind_phi,ind_h) + (chi*(hp - h*(1-d)))^gamma...
                                            - beta*delta*U_hat(ind_phi,ind_hp) - beta*(1-delta)*lambda*opt_p_O(ind_phi,ind_hp)*opt_x_O(ind_phi,ind_hp)...
                                            - beta*(1-delta)*(1-lambda*opt_p_O(ind_phi,ind_hp))*(eta*VO(ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp)); 
                                    end
                                else                                 % N last period - separation shock - does find NOT job from unemployment
                                    simul(id,t).stat = 0;
                                    simul(id,t).wage = 0.0;
                                    simul(id,t).eta = 0.0;
                                end
                        
                            else                   % N last period - no separation shock %
                                % See if successfully searchs OTJ
                                r1 = rand;
                                if (r1 < lambda*opt_p_N(ind_phi,ind_h)) % N last period - no separation shock - successfully searchs OTJ
                            
                                    if (GN(ind_phi,ind_h)==1 || GN(ind_phi,ind_h)==3)             % N - Search OTJ for N
                                        simul(id,t).stat = 1; 
                                        simul(id,t).eta = opt_eta_N(ind_phi,ind_h);
                                        eta = simul(id,t).eta;
                                        ind_hp = opt_ind_hp_N(ind_phi,ind_h);
                                        hp = h_grid(ind_hp);
                                        simul(id,t).wage = eta*(VN(ind_phi,ind_h) - U(ind_phi,ind_h)) + U(ind_phi,ind_h) + (chi*(hp - h*(1-d)))^gamma...
                                            - beta*delta*U_hat(ind_phi,ind_hp) - beta*(1-delta)*lambda*opt_p_N(ind_phi,ind_hp)*opt_x_N(ind_phi,ind_hp)...
                                            - beta*(1-delta)*(1-lambda*opt_p_N(ind_phi,ind_hp))*(eta*VN(ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp)); 
                                    else                                                           % N - Search OTJ for O
                                        simul(id,t).stat = 2; 
                                        simul(id,t).eta = opt_eta_N(ind_phi,ind_h);
                                        eta = simul(id,t).eta;
                                        ind_hp = opt_ind_hp_O(ind_phi,ind_h);
                                        hp = h_grid(ind_hp);
                                        simul(id,t).wage = eta*(VO(ind_phi,ind_h) - U(ind_phi,ind_h)) + U(ind_phi,ind_h) + (chi*(hp - h*(1-d)))^gamma...
                                            - beta*delta*U_hat(ind_phi,ind_hp) - beta*(1-delta)*lambda*opt_p_O(ind_phi,ind_hp)*opt_x_O(ind_phi,ind_hp)...
                                            - beta*(1-delta)*(1-lambda*opt_p_O(ind_phi,ind_hp))*(eta*VO(ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp));    
                                    end
                                else                                         % N last period - no separation shock - does NOT successfully search OTJ 
                                    simul(id,t).stat = simul(id,t-1).stat;
                                    simul(id,t).eta = simul(id,t-1).eta;
                                    eta = simul(id,t).eta;
                                    ind_hp = opt_ind_hp_N(ind_phi,ind_h);
                                    hp = h_grid(ind_hp);
                                    simul(id,t).wage = eta*(VN(ind_phi,ind_h) - U(ind_phi,ind_h)) + U(ind_phi,ind_h) + (chi*(hp - h*(1-d)))^gamma...
                                        - beta*delta*U_hat(ind_phi,ind_hp) - beta*(1-delta)*lambda*opt_p_N(ind_phi,ind_hp)*opt_x_N(ind_phi,ind_hp)...
                                        - beta*(1-delta)*(1-lambda*opt_p_N(ind_phi,ind_hp))*(eta*VN(ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp));
                                end
                            end % End if receives separation shock when N last period     
                        else                                                  % If in O last period  
                            % Recall h from last period and apply any investment
                            ind_h = simul(id,t-1).h_ind; 
                            simul(id,t).h_ind = opt_ind_hp_O(ind_phi,ind_h);
                            ind_h = simul(id,t).h_ind;
                            h = h_grid(ind_h);
                            
                            % If O last period - see if separation shock %
                            rdelta = rand;
                            if (rdelta < delta) % O last period - separation shock    
                                % See if finds job from unemployment
                                r1 = rand;
                                if (r1 < opt_p_U(ind_phi,ind_h))    % N last period - separation shock - finds job from unemployment
                                    if (GU(ind_phi,ind_h)==1 || GU(ind_phi,ind_h)==3)       % Search for N
                                        simul(id,t).stat = 1; 
                                        simul(id,t).eta = opt_eta_U(ind_phi,ind_h);
                                        eta = simul(id,t).eta;
                                        ind_hp = opt_ind_hp_N(ind_phi,ind_h);
                                        hp = h_grid(ind_hp);
                                        simul(id,t).wage = eta*(VN(ind_phi,ind_h) - U(ind_phi,ind_h)) + U(ind_phi,ind_h) + (chi*(hp - h*(1-d)))^gamma...
                                            - beta*delta*U_hat(ind_phi,ind_hp) - beta*(1-delta)*lambda*opt_p_N(ind_phi,ind_hp)*opt_x_N(ind_phi,ind_hp)...
                                            - beta*(1-delta)*(1-lambda*opt_p_N(ind_phi,ind_hp))*(eta*VN(ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp));
                                    else                                                     % Search for O
                                        simul(id,t).stat = 2; 
                                        simul(id,t).eta = opt_eta_U(ind_phi,ind_h);
                                        eta = simul(id,t).eta;
                                        ind_hp = opt_ind_hp_O(ind_phi,ind_h);
                                        hp = h_grid(ind_hp);
                                        simul(id,t).wage = eta*(VO(ind_phi,ind_h) - U(ind_phi,ind_h)) + U(ind_phi,ind_h) + (chi*(hp - h*(1-d)))^gamma...
                                            - beta*delta*U_hat(ind_phi,ind_hp) - beta*(1-delta)*lambda*opt_p_O(ind_phi,ind_hp)*opt_x_O(ind_phi,ind_hp)...
                                            - beta*(1-delta)*(1-lambda*opt_p_O(ind_phi,ind_hp))*(eta*VO(ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp)); 
                                    end
                                
                                else                                 % O last period - separation shock - does find NOT job from unemployment
                                    simul(id,t).stat = 0;
                                    simul(id,t).wage = 0.0;
                                    simul(id,t).eta = 0.0;
                                end
                            else                   % O last period - no separation shock %
                                %See if successfully searchs OTJ
                                r1 = rand;
                                if (r1 < lambda*opt_p_O(ind_phi,ind_h))   % O last period - no separation shock - successfully searchs OTJ
                                    if (GO(ind_phi,ind_h)==1 || GO(ind_phi,ind_h)==3)            % O - Search OTJ for N
                                        simul(id,t).stat = 1; 
                                        simul(id,t).eta = opt_eta_O(ind_phi,ind_h);
                                        eta = simul(id,t).eta;
                                        ind_hp = opt_ind_hp_N(ind_phi,ind_h);
                                        hp = h_grid(ind_hp);
                                        simul(id,t).wage = eta*(VN(ind_phi,ind_h) - U(ind_phi,ind_h)) + U(ind_phi,ind_h) + (chi*(hp - h*(1-d)))^gamma...
                                            - beta*delta*U_hat(ind_phi,ind_hp) - beta*(1-delta)*lambda*opt_p_N(ind_phi,ind_hp)*opt_x_N(ind_phi,ind_hp)...
                                            - beta*(1-delta)*(1-lambda*opt_p_N(ind_phi,ind_hp))*(eta*VN(ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp));  
                                    else                                                          % O - Search OTJ for O
                                        simul(id,t).stat = 2; 
                                        simul(id,t).eta = opt_eta_O(ind_phi,ind_h);
                                        eta = simul(id,t).eta;
                                        ind_hp = opt_ind_hp_O(ind_phi,ind_h);
                                        hp = h_grid(ind_hp);
                                        simul(id,t).wage = eta*(VO(ind_phi,ind_h) - U(ind_phi,ind_h)) + U(ind_phi,ind_h) + (chi*(hp - h*(1-d)))^gamma...
                                            - beta*delta*U_hat(ind_phi,ind_hp) - beta*(1-delta)*lambda*opt_p_O(ind_phi,ind_hp)*opt_x_O(ind_phi,ind_hp)...
                                            - beta*(1-delta)*(1-lambda*opt_p_O(ind_phi,ind_hp))*(eta*VO(ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp));  
                                    end
                                else                                         % O last period - no separation shock - does NOT successfully searchs OTJ
                                    simul(id,t).stat = simul(id,t-1).stat;
                                    simul(id,t).eta = simul(id,t-1).eta;
                                    eta = simul(id,t).eta;
                                    ind_hp = opt_ind_hp_O(ind_phi,ind_h);  
                                    hp = h_grid(ind_hp);
                                    simul(id,t).wage = eta*(VO(ind_phi,ind_h) - U(ind_phi,ind_h)) + U(ind_phi,ind_h) + (chi*(hp - h*(1-d)))^gamma...
                                        - beta*delta*U_hat(ind_phi,ind_hp) - beta*(1-delta)*lambda*opt_p_O(ind_phi,ind_hp)*opt_x_O(ind_phi,ind_hp)...
                                        - beta*(1-delta)*(1-lambda*opt_p_O(ind_phi,ind_hp))*(eta*VO(ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp));  
                                end
                            end % End if receives separation shock when OL last period
                        end % End if for stat (when not dead or died before period start) last period
                    end         % End if dies at end of last period
                end    % End if dead in last period
            end % End if t==1
        end % End loop over time periods in panel
    end % End loop over ids in simul_size
    
    
    count_O =0.0;
    count_N =0.0;
    
    count_E =0.0;
    count_U =0.0;
    for id=1:simul_size
        for t = 1:nt
            if (simul(id,t).stat==2)
                count_O = count_O + 1.0;
            elseif (simul(id,t).stat==1)
                count_N = count_N + 1.0;
            end
            
            if (simul(id,t).stat==0) 
                count_U = count_U + 1.0;
            elseif (simul(id,t).stat==1 || simul(id,t).stat==2)
                count_E = count_E + 1.0;
            end
        end
    end
    
    %% =============================================================================================================== %%
    %========================================== Compute Moments of Interest ============================================%
    %================================================================================================================= %%
    
    % ----------- m1: JTJ Probability
    jtj_rate = 0.0;
    m_sum = 0.0;
    for ind_phi=1:nphi
        for ind_h = 1:nh
    
            ind_hp_N = opt_ind_hp_N(ind_phi,ind_h); 
            ind_hp_O = opt_ind_hp_O(ind_phi,ind_h);
                
            jtj_rate = jtj_rate + lambda*opt_p_N(ind_phi,ind_hp_N)*(1.0-delta)*dist_N(ind_phi,ind_h)*(1.0-nu);
            jtj_rate = jtj_rate + lambda*opt_p_O(ind_phi,ind_hp_O)*(1.0-delta)*dist_O(ind_phi,ind_h)*(1.0-nu);
    
            jtj_rate = jtj_rate + opt_p_U(ind_phi,ind_hp_N)*delta*dist_N(ind_phi,ind_h)*(1.0-nu);
            jtj_rate = jtj_rate + opt_p_U(ind_phi,ind_hp_O)*delta*dist_O(ind_phi,ind_h)*(1.0-nu);
            
            m_sum = m_sum + dist_N(ind_phi,ind_h)*(1.0-nu) + dist_O(ind_phi,ind_h)*(1.0-nu);
        end 
    end 
       
    if (m_sum>0.0) 
        jtj_rate = (jtj_rate/m_sum);
    else
        jtj_rate = 999.0;
    end
    m1 = jtj_rate;

    % ----------- m2: Unemployment Rate
    m2 = (sum(sum(dist_U(:,:)))/( sum(sum(dist_N(:,:))) + sum(sum(dist_O(:,:))) + sum(sum(dist_U(:,:))) ))*100.0;
    
    % ----------- m3: Total Percent Outsourced
    m3 = (( sum(sum(dist_O(:,:))) )/(sum(sum(dist_N(:,:))) + sum(sum(dist_O(:,:)))))*100.0;
    
    % ----------- m4: Ratio Outsourced With to Without Bachelor's Degree
    DO_phi1 = (( sum(dist_O(1,:)) )/( sum(dist_O(1,:)) + sum(dist_N(1,:)) ))*100.0;
    DO_phi2 = (( sum(dist_O(2,:)) )/( sum(dist_O(2,:)) + sum(dist_N(2,:)) ))*100.0;
    m4 = DO_phi1/DO_phi2;
    
    % ----------- m5: Avg. Quarterly Separation Rate
    delta_rate = 0.0;
    m_sum = 0.0;
    
    for ind_phi=1:nphi
        for ind_h = 1:nh
    
            ind_hp_N = opt_ind_hp_N(ind_phi,ind_h); 
            ind_hp_O = opt_ind_hp_O(ind_phi,ind_h);    
    
            delta_rate = delta_rate + (1.0-opt_p_U(ind_phi,ind_hp_N))*delta*dist_N(ind_phi,ind_h)*(1.0-nu);
            delta_rate = delta_rate + (1.0-opt_p_U(ind_phi,ind_hp_O))*delta*dist_O(ind_phi,ind_h)*(1.0-nu);
    
            m_sum = m_sum + dist_N(ind_phi,ind_h)*(1.0-nu) + dist_O(ind_phi,ind_h)*(1.0-nu);
        end 
    end 
    m5 = delta_rate/m_sum;
    
    % ----------- m6: Average Wage Loss After Unemployment Spell   
    wg_uspell = 0;
    m_sum = 0;
    for id = 1:simul_size
        stat_u = 0;
        for t = 1:nt-2 % Need to observe employed wage, unemployment, and then employed wage again
            
            if (simul(id,t).stat~=999)
                if (stat_u==1 && simul(id,t).stat~=0) % Was Unemployed but found Employment
                    w_after = simul(id,t).wage;
                    wg_uspell = wg_uspell + (log(w_after) - log(w_before))*100;
                    m_sum = m_sum +1;
                    stat_u = 0;
                end
    
                if (simul(id,t).stat~=0 && simul(id,t+1).stat==0) % Employed in t and unemployed in t+1
                    w_before = simul(id,t).wage;
                    stat_u = 1;
                end 
            end % End if has not exited the model (if stat~=999)
        end
    end
    wg_uspell = wg_uspell/m_sum;
    m6 = wg_uspell;
    
    % ----------- m7: Wage Dispersion p90/p50
    count_wages = 0;
    for id = 1:simul_size
        for t = 1:nt       
            if (simul(id,t).stat~=0 && simul(id,t).stat~=999)
                count_wages = count_wages + 1;
            end    
        end
    end
    
    % --- First save sample of wages
    wage_sim = zeros(1,count_wages);
    count_wages = 0;
    for id = 1:simul_size
        for t = 1:nt       
            if (simul(id,t).stat~=0 && simul(id,t).stat~=999)
                count_wages = count_wages + 1;
                wage_sim(count_wages) = simul(id,t).wage;
            end    
        end
    end
    
    
    p90_wage = prctile(wage_sim, 90);
    p50_wage = prctile(wage_sim, 50);
      
    m7 = p90_wage/p50_wage;
       
    %----------- m8: Average Annual Real Wage Growth
    % average wage growth over 1 year
    wg_sum2 = 0.0;
    m_sum = 0.0;
    for id = 1:simul_size
        for t = 1:(nt-4)
            if (simul(id,t).stat~=0 && simul(id,t+4).stat~=0 && simul(id,t).stat~=999 && simul(id,t+4).stat~=999) 
                w1 = simul(id,t).wage;
                w5 = simul(id,t+4).wage; 
                wg_sum2 = wg_sum2 + (log(w5) - log(w1))*100.0;
                m_sum = m_sum + 1.0;
            end 
        end 
    end 
    wg_sum2 = wg_sum2/m_sum;
    m8 = wg_sum2;
    
    % ----------- m9: College Wage Premium
    avg_wage_phi1 = 0.0;
    avg_wage_phi2 = 0.0;
    sum_phi1= 0.0;
    sum_phi2= 0.0;
    
    for id = 1:simul_size
        for t = 21:nt % Age >= 25 yrs
            
            if (simul(id,t).stat~=0 && simul(id,t).stat~=999)   % If Not Unemployed or dead
                if (simul(id,t).phi_ind == 1) 
                    avg_wage_phi1 = avg_wage_phi1 + simul(id,t).wage;
                    sum_phi1 = sum_phi1 + 1.0;
                else
                    avg_wage_phi2 = avg_wage_phi2 + simul(id,t).wage;
                    sum_phi2 = sum_phi2 + 1.0;
                end 
            end 
        end 
    end 
    
    avg_wage_phi1 = avg_wage_phi1/sum_phi1;
    avg_wage_phi2 = avg_wage_phi2/sum_phi2;
    
    m9 = avg_wage_phi2/avg_wage_phi1;
    
    
    %---------- m10 - Unemployment Rate Among Age 20
    count_E =0.0;
    count_U =0.0;
    for id = 1:simul_size
        t=1; % age 20
        if (simul(id,t).stat~=999)
            if (simul(id,t).stat==0)  
                count_U = count_U + 1.0;
            elseif (simul(id,t).stat==1 || simul(id,t).stat==2) 
                count_E = count_E + 1.0;
            end    
        end 
    end 
    m10 = (count_U/(count_E + count_U))*100.0;
    
    mm_D1 =[m1;m2;m3;m4;m5;m6;m7;m8;m9;m10]; % Model Moments

end % end function